Bose— Einstein condensation and the magnetically ordered state of TICUCI3 
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The dimerized S = \ spins of the Cu 2+ ions in TICUCI3 are ordered antiferromagnetically in the 
presence of a field larger than about 54 kOe in the zero-temperature limit. Within the mean-field 
approximation all thermal effects are frozen out below 6 K. Nevertheless, experiments show signif- 
icant changes of the critical field and the magnetization below this temperature, which reflect the 
presence of low-energetic dimer-spin excitations. We calculate the dimer-spin correlation functions 
within a self-consistent random-phase approximation, using as input the effective exchange coupling 
parameters obtained from the measured excitation spectra. The calculated critical field and magne- 
tization curves exhibit the main features of those measured experimentally, but differ in important 
respects from the predictions of simplified boson models. 

PACS numbers: 75.10.-b, 75.30.-m, 67.85.Jk 
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I. INTRODUCTION 

The concept of Bose-Einstein condensation dates back 
more than 80 years to the prediction of Einstein, based 
on Bose's work on the statistics of photons, that a gas 
of non-interacting massive bosons would condense below 
a certain critical temperature T c . The condensation im- 
plies that below T c a non-zero fraction of the total num- 
ber of particles occupies the lowest single- particle quan- 
tum state. For dilute atomic gases this phenomenon was 
realized experimentally in 1995 for trapped clouds of al- 
kali atoms (see e.g. Ref. [l|). 

For a uniform gas of density n the transition temper- 
ature is given by kT c « 3.31?i 2 n 2 / 3 /m, where m is the 
particle mass. For particles trapped in a harmonic oscil- 
lator potential (trap frequencies u> x , uj y and oj z ) one has 
kT c w 0.9471(7^^1^ ) 1/3 , where N is the total number 
of particles. In the latter case the particle mass enters 
through the trap frequencies, equal to the square root of 
the force constants in the three directions divided by the 
particle mass. When a trapped gas is dilute in the sense 
that the atom-atom scattering length is much less than 
the interatomic distance, the observed transition temper- 
atures agree well with theoretical expectation for a non- 
interacting gas. For less dilute gases interaction effects 
give rise to an observable small shift of T c proportional 
to the scattering length. 

The condensation of massive bosons into a single quan- 
tum state is intimately connected to the conservation of 
particle number. For massless bosons such as phonons 
or magnons the particle number is not fixed but depends 
on temperature, and there is therefore no Bose-Einstein 
condensation in the traditional sense of the term. How- 
ever, there has been a wide use of model Hamiltoni- 
ans for magnetic systems that have features in common 
with those of interacting, massive bosons. The aim of 
the present work is to consider one such specific system, 
that of the dimerized Cu 2+ spins in TICUCI3, and com- 
pare predictions of such models with calculations that are 
based on (approximate) solutions of the many-body prob- 
lem of interacting spins. A recent review of experimental 



and theoretical developments concerning the magnetic 
ordering of TICUCI3 and related compounds has been 
given by Giamarchi et al£ 

Magnetization measurements^ and inelastic neutron- 
scattering experiments^ demonstrate clearly that 
nearest-neighboring pairs of S = i spins of the Cu 2+ 
ions in TICUCI3 are dimerized leading to an S — ground 
state and an S — 1 excited triplet around 5.2-5.7 meV 
above the singlet. Due to the exchange interactions be- 
tween the dimers the excitations become strongly disper- 
sive, and, in the zero-temperature limit, the minimum 
energy of the degenerate singlet-triplet mode, at (001), 
is only about 0.7 meV. When a field is applied, the energy 
of one of the three normal modes is reduced and goes to 
zero at a critical field of about 54 kOe at zero tempera- 
ture. The Cu spins of isostructural KCUCI3 are similarly 
dimerized, but the interdimer interactions are relatively 
weaker and the critical field is about 230 kOe at T = in 
this system* 7 " The phase transition shown by TICUCI3 at 
the field where the excitation energy vanishes, has been 
analyzed by Nikuni et alM- They assumed the dimer sys- 
tem to be described by an effective Haniltonian of the 
form 
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where the bosonic operators a* and a denote "magnon" 
creation and annihilation operators and the positive con- 
stant vq denotes the strength of the repulsive magnon- 
magnon interaction, assumed to be a delta function in 
real space. The quantity /1 plays the role of a chemical 
potential, assumed proportional to the difference between 
the applied magnetic field and the critical field. Using vq 
and m as fitting parameters, Nikuni et al. were able to 
give a reasonable account of the temperature dependence 
of the critical field and the magnetization along the ap- 
plied field in the ordered antiferromagnetic state. An 
extended version of their theory based on a more realis- 
tic dispersion of the magnetic excitations was presented 
by Misguich and Oshikawa. 9 However, as we shall see in 
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detail in Sec. IV below, these simplified boson models 
suffer from inconsistencies that originate in their neglect 
of the highest level of the triplet. Another, more general 
problem with boson models is that double occupancy of 
a local site should be prohibited. In this connection we 
mention the work of Sirker et a/.,— who used a bond- 
operator approach to map the spin system onto a model 
of interacting bosons by introducing an infinite on-site 
repulsion between local triplet excitations. 

In the following treatment of the dimcrizcd spin system 
in TICUCI3 we adopt a different point of view and start 
from the Hamiltonian for the spin system itself, using as 
input the effective exchange coupling that has been de- 
rived from measured excitation spectra. Our approach is 
a generalization of the zero-temperature theory by Mat- 
sumoto et al.^ who used the random-phase approxima- 
tion (RPA) to calculate magnetization curves and exci- 
tation spectra for TICuCL;. They also considered the 
case when the phase transition is induced by the applica- 
tion of a hydrostatic pressure. 12 ^ 3 Here we only address 
the case of a field-induced transition, but the theory of 
Matsumoto et al. is extended to include both the effects 
of quantum fluctuations at zero temperature and the ef- 
fects of thermal fluctuations. The self-consistent version 
of the RPA, which is the one applied here, is faced with 
the similar problem of double occupancy as the boson 
modelling. However, here this problem is found to have 
a natural solution by a consideration of the higher order 
modifications of the Green functions. The self-consistent 
RPA theory for the paramagnetic phase of the dimer-spin 
system is presented in Sec. II, which, in Sec. Ill, is fol- 
lowed by an analysis of the antiferromagnetic phase. A 
closer examination and discussion of the results obtained 
are referred to the last Sec. IV. 



II. EXCITATIONS IN THE PARAMAGNETIC 
PHASE 

A. The self-consistent RPA theory 

The TICuCL crystal is monoclinic (space group P2\jc) 
and the lattice parameters are a — 3.9815 A, b = 14.144 
A, c = 8.8904 A and /3 = 96.32° at room temperature^ 
The crystal is constructed from layers with configuration 
CU2CI6 stacked on top of each other so as to form two 
chains of Cu ions parallel to the a axis. The chains are 
separated by Tl ions and pass through the center and 
corners of the b-c plane in the unit cell. There are four 
Cu ions or two dimer pairs per unit cell. The dimer 
pair in the unit cell belonging to the chain through a 
corner is located at site 1: (x,y,z) and site 2: (x, y, z), 
and the pair belonging to the other chain is placed at 
site 3: (x,y + \,z + |) and site 4: (x,y + |,z + ^). 
Here x = 0.2338, y = 0.0486, and z = -0.0175, and 
the numbering of the sites from 1 to 4 defines the four 
different Cu-sublattices. 



The Hamiltonian is assumed to be 

M^i^^.^-^^if.^, (2) 

ij i 

where s i is the spin-variable of the Cu ion at the ith 
site. The most important exchange parameter is A = 
where (1^2) are the nearest-neighbor Cu pairs 
(the 1-2 or the 3-4 ions in the unit cell). The Fourier 
transform of the Heisenberg exchange interactions be- 
tween spins on sublattice a and (3 is defined in terms of 
the remaining coupling parameters 

M«)= £' m)z- iq<Ri - Ri \ ieo-suU., (3) 

jS/3-subl. 

where is the position of the ith dimer, and the prime 
indicates that the dominating interaction A is excluded 
from the sum, (ij) ^ (iiia)- 

When the interactions between the dimers are ne- 
glected, the Hamiltonian may be diagonalized exactly 
in terms of independent products of single-dimer eigen- 
states. The total spin of the zth dimer is S i = + Sj 2 , 
where and Si 2 denote the Cu spins belonging to, re- 
spectively, the sublattices 1 and 2, or 3 and 4. The total 
spin defines the basis \SS Z ), and when the field is along 
the z axis, the eigenstates of the non-interacting dimer 
are 

state |3) = |1 — 1) at the energy A + h, 

state |2) = 1 10) at the energy A, 

state |1) = |1 + 1) at the energy A — h, 

state |0) = 1 00) at zero energy, 

with h = g/j, B \H\ and A positive. When h < A the 
ground state is the non-magnetic singlet |00). In the 
present section we focus on this condition, and we shall 
assume that the system stays paramagnetic also in the 
presence of the interdimer interaction a «(g). The orig- 
inal Hamiltonian |(5J| may be rewritten in terms of two 
dimer-spin variables, the sum Si = + Si 2 and the dif- 
ference Si = Sjj — Si 2 . The sum operator only has non- 
zero matrix elements between the three excited S = 1 
states. We consider the case where the populations of 
these levels are small (at sufficiently low temperatures in 
the disordered phase), in which case the dynamical ef- 
fects due to Si are negligible. When Si is neglected, the 
Hamiltonian involves only a single effective q-dependent 
exchange term — -| J2 q J(q)Sq ' S-q with 

J(q) = \ [3u(q)+ 322(g) -2u(q)- 021(g)] 

± \ \3xz{q) + 024 (g) - du(q) - 023 (<?)] • (4) 

The presence of two equivalent dimers per unit cell yields 
two values for the effective interaction for each value of 
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the wave vector within the first Brillouin zone. Alter- 
natively, one may use an extended zone scheme with an 
effective basis of one dimer per unit cell, in which case 
only the upper sign applies. 

In order to study the spin dynamics of this Hamilto- 
nian we introduce the standard basis operator s] 15 ' 16 ! 17 for 
the jth dimer 



a U = (lM>H)j > M» v = 0, 1, 2, 3. 



(5) 



In the present case of a dimer system with stationary 
bonds, these operators serve the same purpose but are 
of more general use than the "bond operators" applied 
by Matsumoto et al. 11 In terms of the standard basis 
operators the components of Sj become 
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and the Hamiltonian may be written 
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when only the Fourier transform J(ij) of the effective 
interaction, Eq. is included. Next we define a 6 x 6 
matrix of Green functions^ 



G(ij,w) 



0(t)([aM,aUo)])e^dt. (8) 



A single bracket (• • • ) denotes the thermal expectation 
value, and a^t) is a vector operator of site i at time t 
with components 



a 01; a lCb a 02; a 20; a 03> a 3oJ ■ L y J 

With the short-hand double-bracket notation G(ij, u>) — 
((djjat)), the equations of motion for the Green func- 
tions are (see for instance Ref. [l7t ) 

fkj((a i; a])) - (([a,, %} ; a])) = ( [a t , a}] ) , (10) 

where the new higher-order Green functions introduced 
by the second term are determined by the Hamiltonian 
([7|) by the use of the commutator relation 

a i>„>] = S ij (<V a ^' - V' a Mv) • (11) 

By using an RPA decoupling of the higher-order Green 
functions, a^a 3 ^ v , » a^(a^ v ,) + (a^,)a£ v , one finds 
that the equations of motion are reduced to a closed set 



of equations, which are solvable by a Fourier transforma- 
tion. Further, the 6x6 matrix equations decouple into 
3 sets of 2 x 2 matrix equations, and one of these is 



Tiw-E, -J, \ (G n G w 
J 3 hu + E 3 ) \G 61 G 66 



n m 
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03 



. (12) 

The Green functions depend on the Fourier variables, 
Gn V — G„ u (q,cj), and we have introduced the following 
parameters 

E 1 = A-h-J 1 , J 1 = n 01 J(q), 

E 3 = A + h-J 3 , J 3 = n 03 J(q) , (13) 

with being the average population of the /ith dimer 
level 



n„ = (a!j„) , and the difference n^v 



%~n u . (14) 
Inversion of the first matrix in Eq. (|12p results in 
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n m {E 3 + hw) 



n 03 J 1 
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The poles of the Green functions determine the excitation 
energies E q , which are given by 



E± =E q ±(h- \n 13 J{q)) , 



(16) 



E 2 q = A 2 ~ (n 01 + n 



,)AJ(q) + (in 13 J(q)) : 



The (2,5) part of G(q, w) is given by the same expression 
as the (1,6) part in Eq. (fl~5|) . except that u> is replaced by 
— oj. The result for the (3,4) part is that obtained from 
Eq. (|15p for h = 0, corresponding to the replacement of 
n 01 and n 03 by n Q2 , and it leads to poles at the energies 



±E*, where 



E\ = VA 2 - 2A J 2 , J 2 = n 02 J(q) 



(17) 



We introduce the following matrix of equal-time corre- 
lation functions 



I( q ) = i^(a i a]>e-^-^ 



N 



(18) 



where N is the number of dimers. According to the 
fluctuation-dissipation theorem (see e.g. Ref. [13) 



A{q) 
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G"(q,u)d(hu)), (19) 



where G" denotes the imaginary part of the matrix Green 
function and (3 = 1/kT. By definition, the average value 
of, for instance, the 11-component is 



ooi °lo) = (aoo) = »o ■ ( 2 °) 



x 00/ 
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Hence, by calculating the q averages of the correlation 
functions, the Green functions may be used for deter- 
mining the populations of the dimer levels. From the 
diagonal part of the A-matrix, we get straightforwardly 



n Q + m ^n + n 3 
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(22) 



(23) 



The three equations determine the four population num- 
bers, when they are supplemented by the exact condition 
that 



•7i3 = l 



(24) 



The RPA decoupling is valid in the mean-field (MF) 
approximation, where the thermal averages are de- 
termined by the Hamiltonian !K MF for the non- 
interacting system, i.e. within the approximation {a^J) 
Tr[a My exp(— /3J{ MF )] = 5^ v n^ F . In the case of strong 
dispersion, this approximation certainly underestimates 
the populations of the excited levels. Anticipating that 
the correlation effects predicted by the RPA theory 
are reasonably trustworthy, the self-consistent equations 
above should lead to a more accurate determination of 
the population numbers than that offered by the MF 
approximation. Unfortunately, the present theory also 
predicts that the off-diagonal components of the A ma- 
trix are non-zero contradicting that, for instance, A 16 = 
(a l 01 a Q3 ) should vanish identically. Phrased differently, 
this inconsistency implies that the RPA result for the 
occupation numbers is not unique but depends on the 
particular choice of correlation functions used in the cal- 
culation. A non-zero value of (a ia 03 ) is the equivalent 
of a double occupancy of bosons at a single site. In- 
stead of introducing an arbitrary repulsive potential, we 
are here going to consider possible improvements of the 
RPA-decoupling procedure applied above. 

The present system is in many ways similar to the 
singlet-doublet system encountered in praseodymium 
metal. Improvements of the RPA for this system have 
been derived in Ref. [l^, and this theory has been applied 
to the calculation of the linewidths and the energy renor- 
malization of the excitations in two other dimer systems 
Cs 2 Cr 2 Br9 and KCuCl 3 for the case of zero field It is 
straightforward to generalize the singlet-doublet theory 
so as to account for the presence of a third level. Here, we 
are going to extend the theory to the case where the field 



is non-zero, and, in the next section, to consider the mod- 
ifications produced by an ordered moment. The presence 
of the intcrdimer interactions J(ij) implies that Og does 
not commute with the Hamiltonian. Hence, the assumed 
ground state, the product state of |0)i, is not an eigen- 
state of the interacting system. The situation compares 
with the simple Heisenberg antiferromagnet, where the 
mean-field Neel state is not the true ground state. In this 
case the RPA theory predicts a zero-temperature reduc- 
tion of the antiferromagnetic moment from its saturated 
Neel state value. Equivalently, the RPA results above 
imply that n is smaller than its saturation value 1 at 
zero temperature. The single-dimer population numbers 
are subject to quantum as well as thermal fluctuations. 

One of the terms neglected in the RPA equation (jT3J) 
involves the Green function (((a 00 — a l n — n 01 )a 01 ;&l.}). 
Since a 00 and are not true constants of motion, this 
Green function may modify the RPA result. The equa- 
tions of motion for the Green functions neglected in Eq. 
(fT2| have been analyzed in Ref. [l^. The consequences are 
that the RPA parameters are being replaced by effective 
ones and that the excitations become damped. One of 
the effective parameters J xy (q) replaces J(q) in J 1 and 
J 3 in all of the equations above and 



Jxyil) =J(*l)- a xy + Vxy [hyi") ~ O-xy] 



where 



n. 



1 . 



(25) 



(26) 



The most important modification of J(q) is the constant 
shift introduced by a xyl but first we want to discuss the 
other term 



hy{w) 



(27) 



where Xxyi^^) ^ s a generalized susceptibility^ The 
imaginary part of b xy (uj) determines the damping effects, 
which are, however, small at low temperatures. The 
renormalization of J xy (q) produced by the real part of 
b xy {u>) is somewhat smaller than the constant shift due 
to a xy , but it is not entirely negligible. Because of its 
moderate importance we have simplified the expression 
for the b xy (ui) term as follows 



b xy (u>) » Re [b xy (E q /h)] « (j 



'01 



B n 



[J(fe)] 2 [J(q)-J(fc)] 



[J(q) 



J(k)f 



(28) 



We include only the real part and neglect modifications 
produced by the field. The quantity e is due to the finite 
lifetime of the excitations and we take it to be a con- 
stant, e w 0.1 meV at zero temperature and field. The 
linewidths are going to increase rapidly when the temper- 
ature becomes comparable to A/fc. We have accounted 
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for this effect in a rough manner by scaling the result 
B q at zero temperature and field by the population- 
dependent factor in front, where n° is the value of a pop- 
ulation number at zero field and temperature. A closer 
examination indicates that this simple scaling accounts 
for the increase of the linewidths in a reasonable way. 
The scale factor is unimportant for the analysis of the 
low temperature properties, but the power of 3 used in 
this expression ensures that the scale factor times r] xy 
vanishes in the high-temperature limit. 

Finally, the most important renormalization effect, the 
constant term (f + r\ xy )a xy in J xy (q) in Eq. (j2"B"|) , is de- 
termined implicitly by 

ly J xy (q) ( 1 + n - +n +) =Q (29) 

N q E q + E q 

The renormalized value of A 16 is equal to this sum over q 
times n 01 n 03 , and since the sum now vanishes, the condi- 
tion (aQiCLQ 3 ) = is satisfied. The parameters J 1 and J 3 
in Eqs. (flli)) - (fTrJ)) are replaced by, respectively, n 01 J xy (q) 
and n 03 J xy (q), and, similarly, J 2 in Eq. (|17p is replaced 
by n 02 J z (q), where the effective exchange coupling is 

J z (q) = J{q)-a z + Vz [b z {u>)-a z ] ■ (30) 

Here rj z = l/ng 2 - 1 and b z (u) w (n 02 /n^ 2 ) 3 B q . The 
constant term (f + rj z )a z is determined by the condition 
that (<2o2 a o2) = 0. Besides the modifications of the ex- 
change couplings, the energy splitting A in E 1 and E 3 , 
or in E 2 , is replaced by, respectively^ 

A 1 = A 3 = A + (o iB1 , + o,)/2, 

A 2 = A + a xy . (31) 

The renormalization of the energy-level separations has 
not much influence on the final calculations, and the 
small additional modifications derived in Ref. might 
have been neglected. However, the leading order effects 
of these extra terms are actually included above by as- 
suming 1 + r/ xy to be a factor n + (n 1 + n 3 )/2 smaller 
than derived in Ref. fl9L and, similarly, 1 + i] z has been 
divided by n + n 2 . 

B. Comparison with experiments 

Cavadini et al. 5 and Oosawa et al& have measured the 
dispersion of the magnetic excitations in TICUCI3 at zero 
field in the zero-temperature limit (1.5 K). The two sets 
of results agree where they overlap, and the combined 
experimental results are most closely reproduced by the 
dispersion parameters derived by Oosawa et al. The pa- 
rameters used here (in units of meV) determine the ef- 
fective exchange coupling according to 

J eS (q) = 0.46cos(g-a)-0.05cos(2g-a) 
+ 1.53 cos (q ■ (2a + c)) 
T 0.86cos(q-(a + ic))cos(iq-6) . (32) 
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FIG. 1: (Color online) The susceptibility of TlCuCl 3 deter- 
mined experimentally with a field of 10 kOe along the 6 axis. 
The pluses show the data of Oosawa et a/.,— and the open 
circles are the results of Dell'Amore et al£ The experimental 
results are compared with the calculated ones obtained by as- 
suming g = 1.97 or g — 2.33. In all other calculations we use 
g = 2.06. 

These are the parameters derived by Oosawa et al. except 
that their coupling between the two chains 
=F[0.98 cos(q • (a + |c)) - 0.12 cos(iq • c)] cos(iq • b) has 
been approximated by a single term. The three modes 
are degenerate at zero field, and the dispersion relation 
assumed by Oosawa et al. in their analysis is E q = A^ ff — 
2A cff J off(<?) implying that 

A 'ff = E E l = ( A + a o) ( A + °o + 2a Ki) . 
q 

J cS (9) = < [J(q) + 4 v B q ] . (33) 

Here a is the value of a xy or a z at zero temperature 
and field. The quantity B q is roughly proportional to 
J(q) and its averaged value with respect to q is zero. 
The RPA equations above have been solved numerically 
by an iterative procedure. The calculations benefit from 
the fact that all q summations may be parameterized 
in terms of J eff (q) [we neglect the minor difference be- 
tween J(q) + VxyB q and J(q) in Eq. (|2"5)) determining 
B q \. Hence all summations may be expressed as inte- 
grals with respect to J eS (q) times a corresponding "den- 
sity of states" calculated once and for all from Eq. (|32[) . 
The value of A eff used in the calculations is 5.671 meV, 
(almost) equal to the value of 5.68 meV derived by Oo- 
sawa et al. 6 Besides this parameter and those defining 
J off (g) we have assumed that g = 2.06, which is the gen- 
erally accepted value for g in the case where the field 
is applied along the b axisi^i22 Finally, we have added 
the mean field from the parallel component to the ap- 
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FIG. 2: The solid line is the theoretical result for the critical 
field as a function of temperature using the exchange param- 
eters introduced by Eq. (|32|l . The dashed line is the result 
obtained if using instead the exchange parameters of Oosawa 
et al£ The critical field is here defined to be the one at which 
the paramagnetic phase becomes unstable. The experimental 
points are those obtained when the field is applied in the b 
direction by Oosawa et al. 2 ^ and Shindo and Tanaka^ 
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3: (Color online) The parallel magnetic moment per Cu 
ion, m z , as a function of temperature calculated at various 
values of the field applied along the b direction. In the case 
of H — 53 kOe the system is predicted to stay disordered all 
the way to zero temperature. The experimental points are a 
selection of those obtained by Oosawa et a/.— £ 



plied field so that gfi B H = h in Eq. |(3J) is replaced by 
h = h + J F (0)(S z ). The ferromagnetic coupling J F (0) 
is estimated to be about —1.9 meV by Dell'Amore et al~ 
and of the order of —2.8 meV by Oosawa et alM (using 
their parameters determined by a cluster series expan- 
sion). Here we assume Jp(0) = —2.4 meV. The moment 
per Cu 2+ ion parallel to the applied field is 

1 \ - 1 n X — Tin 

m * = ^Z^BTJ^/ = Mb 2 ' ( ' 

i 

It is worthwhile to notice that although the population 
numbers of the excited levels are predicted to be non- 
zero at T = 0, Eq. (j2"T]) - (f2"4")) . the quantum fluctuations 
do not give rise to any difference between n 1 and n 3 , i.e. 
m z — at zero temperature as long as the system stays 
paramagnetic. This is consistent with the condition that 
Si ^iz = ( a i i — a 33) commutes with the Hamiltonian. 

Using the model defined here we have calculated the 
susceptibility as a function of temperature. The result 
is compared with experiments in Fig. [TJ The calculated 
critical field at which the paramagnetic phase becomes 
unstable is compared with experiments in Fig. [3J The 
induced magnetic moment m z at various values of the 
field has been calculated as a function of T. These results 
are shown in Fig. [3) This figure also includes results 
obtained in the ordered phase, which is considered in the 
following section. 



III. EXCITATIONS IN THE 
ANTIFERROMAGNETIC PHASE 

The paramagnetic phase becomes unstable when the 
energy of the lowest excitation vanishes. The lowest en- 
ergy mode is the one with the energy E~ at q = Q = 
(001), and below the transition the expectation value of 
the dimer-spin variable S i becomes non-zero. The or- 
dering is antifcrromagnetic in the sense that (S^ have 
opposite signs on the two chain sublattices. The antifer- 
romagnetic ordering may be transformed to the uniform 
one by an interchange of the two spins in the definition 
of S i for the dimers belonging to, for instance, the 3-4 
sublattices. The only effect of this transformation is that 
the interchain coupling between the 1-2 and 3-4 dimers 
changes sign, i.e. the ± in front of the second term of 
J(q) in Eq. (01 is being replaced by =F, and J(Q) and 
J(0) are being interchanged within the extended zone 
scheme. That the ordering is antiferromagnetic instead 
of being uniform does not introduce any further compli- 
cations. (S^ is perpendicular to the z-direction of the 
field, but its direction within the x-y plane is arbitrary 
(as long as any anisotropy is neglected). For convenience 
we shall make the choice that the ordered moment is 
along the x axis, and we define 

m xy = gn B ^, (S x ) = lj2{S ix )e^ . (35) 

i 

When (S x ) is non-zero, the MF Hamiltonian for the 
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"non-interacting" ith dimer becomes 



= (A - h) a\ ± + A a\ 2 + (A + h) a 33 
1 

'V2 



J(Q)-= [4 - a\ + a* 3 - ai x ] (S x ) (36) 



in terms of the standard basis operators of the paramag- 
netic system. We shall continue to label the eigenstates 
by |/x), where fi = 0, 1, 2, 3, and the ground state of this 
MF Hamiltonian may then be written 

|0> = cos(9|00) -sin 6» [cos (a - f) |1 + 1) 

+ sin(a-f) |1-1)] . (37) 

The two angles 9 and a minimize (0|:Hm F |0), or may be 
determined by demanding the off-diagonal terms of the 
MF Hamiltonian to vanish, and we get 



tana = 



h cos 29 



and 



sin26> = 



A cos 2 6 — h sin 2 9 sin 2a 
2J(Q) cos a 



(38) 



A cos 2a 



(cos 26» + 2 sin 2 8 sin 2 a) (S x ) 



(39) 

By calculating the state vectors of the excited MF levels 
to order 9 2 , we find that the order parameter is 

(&x) = ~~ ^ sm 2^cosa — n 13 sin6'sina + 0(9 3 ) , 

where the higher-order terms 0(6 |3 ) vanish if n = 1, and, 
likewise, the ferromagnetic component is 



<"13 ' 



"oi + n Q3 . 2 . nra^ tA-w 
sin sin 2a + (j (9 ). (41) 



By h c we denote the critical field at which these equations 
have a non-zero solution for (S x ) in the limit of 9 — ► 0, 
and this field is found to be determined by 



1 - 



A 



J(Q) 



n 01 + n 03 



'13 



(42) 



If we replace J(Q) by J xy {Q) and A by A x this condition 
is the same as that derived from the requirement that the 
energy of the lowest paramagnetic excitation Eq , within 
the self-consistent RPA, should vanish at the transition. 
The results above are more general but coincide with 
those derived in the zero-temperature mean-field theory 
of Matsumoto et al 11 [Their angle <f> corresponds to our 
j — a. In our notation h = h Q + J F (0)(S z ), whereas in 
their notation h in Eq. (|38|) should read h — J(Q)(S Z ) 
corresponding to the replacement of h by h with their 
implicit assumption that the ferromagnetic interaction 
J F (0) is equal to —J(Q)}. 

When the MF Hamiltonian of the antifcrromagnetic 
phase has been diagonalized we may proceed as in the 
paramagnetic case for calculating the correlation func- 
tions. The positions of the four different levels and the 



matrix elements of S i may be calculated analytically if 
only terms to leading order in 9 2 are included. We are 
not going to present these results, since in the final cal- 
culations we chose the more accurate approach of diag- 
onalizing the MF Hamiltonian numerically. In terms of 
the standard basis operators of the final MF Hamilto- 
nian, Sj X is now [m 3 (a3 + a J 03 ) ~ m 1 (a{ () + Ooi)]/v2) 
where m 1 and m 3 are different from 1 and from each 
other. Here we neglect the extra complication that the 
matrix elements between the excited states, as for in- 
stance (l|5ja;|3) oc 9, become non-zero. These excited 
state contributions to the correlation functions get mul- 
tiplied by n 13 , hence they are unimportant not only when 
the ordered moment is small, but in most of the regime 
where the RPA modifications of the MF behavior are of 
importance. 

The final Hamiltonian may be written in the same way 
as in the disordered case, Eq. ([7]) . The positions of the 
three excited levels are being shifted and J(ij) is being 
multiplied by different factors depending on which op- 
erator product is considered, and, finally, the remaining 
off-diagonal products, ctoi a oi' a oi a 30 e ^ c -' now apP ear - 
The new off-diagonal contributions are all of order 9 2 and 
affect the diagonal correlation functions only to order 9 A . 
Hence, to leading order these extra contributions may be 
neglected. In this case the matrix equations once again 
decouple into 3 sets of 2 x 2 equations, which may be 
solved analytically. The result for the population num- 
bers is the same as that given by Eqs. (|2"l"T) - (f2"4")) . except 
that A and A ± h are being replaced by the energies of 
the three corresponding MF levels and that J t , J 2 , and 
J 3 are being multiplied by matrix-element factors, which 
are slightly different from 1 and from each other. The 
equivalence implies that the modifications of the RPA 
correlation functions may be calculated as in the para- 
magnetic case, and, for instance, a xy is still determined 
by Eq. (|2"9"| except that the expressions for are be- 
ing modified. We used this approximation, valid in the 
limit of 9 2 being small, for calculating the magnetization 
curves. The results were close to those shown in Fig. |31 
however, in the final calculations we included the higher- 
order modifications. 

For a given set of population numbers the MF Hamilto- 
nian was diagonalized numerically, determining all pos- 
sible matrix elements and the four energy levels. The 
knowledge of the population numbers and the matrix el- 
ements is used for determining the two expectation values 
(S x ) and (S z ), and for constructing the total Hamiltonian 
expressed in terms of the standard basis operators. When 
the interactions between the excited states are neglected, 
the equations of motion lead to a 4 x 4 set of matrix equa- 
tions for the xy part and a 2 x 2 set for the longitudinal 
Green functions. The two sets of equations were inverted 
analytically (utilizing Mathematica for handling the set 
of 4 x 4 equations). In this way we derived an explicit 
expression for the correlation function matrix A(q). The 
averaged values of the diagonal components were used for 
calculating the population numbers as in the paramag- 
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5 




Magnetic field (kOe) 

FIG. 4: (Color online) The minimum energies of the three 
different dimer excitations as functions of field at 1.5 K. The 
calculated results (the solid lines) are compared with the ex- 
perimental results of Riiegg et aZ . 25 ' 26 



netic case. The renormalization parameter a xy and a z are 
determined by, respectively, A 16 = and A 34 = 0. The 
condition A 16 = also implies that A 25 — 0, but not nec- 
essarily that the new off-diagonal components vanish. We 
have neglected the possibility that the renormalization of 
the additional off-diagonal exchange terms might be dif- 
ferent, since the new terms are small whenever the renor- 
malization effects are important. Except for the matrix- 
element modification of the exchange terms we use the 
same approximate expression, Eq. (|28[) . for the renormal- 
ization parameter b (u>), and similarly for b z (u>). When 
the renormalization parameters have been determined we 
may calculate the renormalized value of J(Q) in the MF 
Hamiltonian (j36|) , which is being replaced by J xy (Q) de- 
termined from Eq. Similarly A ± h and A in Eq. 
([Ml) are replaced by, respectively, A 1 ±h and A 2 given by 
Eq. J2IJ). The MF Hamiltonian is then consistent with 
the renormalized RPA expressions, and the whole proce- 
dure has been carried out in a self-consistent manner, so 
that the population numbers assumed as a start are the 
same as those derived. 

This theory was used for calculating the properties of 
the dimer system in the ordered phase. The temperature 
dependence of the parallel moment in a constant applied 
field is shown in Fig. [3] fn Fig. Q] we show the calculated 
energies for the three modes at (001) as functions of field 
at 1.5 K compared with the neutron-scattering results 
of Riiegg et a/. 25 ' 26 The experimental results indicate 
that A for the longitudinal mode, the energy of which is 
nearly unaffected by the field in the paramagnetic phase, 
is slightly larger than for the two other modes, and we 
have accounted for this effect by adding 0.03 meV to 
A 2 in Eq. ([3T]) . This modification does not affect the xy- 




Magnetic field (kOe) 

FIG. 5: (Color online) The squares are the experimental re- 
sults for the field dependence of the parallel moment m z at 
1.3 K and are the data of Tatani et al^ presented in Ref. 
ITU . The open circles are the neutron diffraction results for 
the square of the ordered antiferromagnetic moment m% y ob- 
tained by Tanaka et alM- with the field along the b direction 
at 0.2 K. The solid lines show the corresponding theoretical 
predictions, and the dashed ones are the results of using the 
MF approximation. 

polarized modes, and as discussed by Matsumoto et aL,— 
the lowest-energy xy-polarized mode becomes the Gold- 
stone mode in the ordered phase, the energy of which de- 
pends linearly on | q — Q\ and is zero at the ordering wave 
vector^ The spin-resonance experiments of Glazkov et 
al£^ indicate that the Goldstone mode develops an en- 
ergy gap for fields larger than the critical one correspond- 
ing to the presence of a small anisotropy within the x—y 
plane of the same order of magnitude as the anisotropy 
considered above. The two anisotropy terms are unim- 
portant for the renormalization effects and are neglected 
elsewhere in the present calculations. The final figure, 
Fig. shows the field dependencies of the squared pri- 
mary order parameter m xy and of the parallel magnetiza- 
tion m z in the zero-temperature limit. The comparison 
with experiments shows that the self-consistent RPA ac- 
counts reasonably well for the field dependence of m z , 
whereas the order parameter is calculated to increase 
rather faster with field than observed. The corresponding 
MF model, on the other hand, underestimates the value 
of m z , but predicts an order parameter which is close to 
the one observed. 



IV. DISCUSSION 

The present dimer system is unique because it clearly 
exhibits the importance of quantum fluctuations. The 
system is close to a quantum critical point, and a zero- 
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FIG. 6: (Color online) The renormalization parameters a xy 
and a z (left scale) and the population numbers (right scale) 
calculated as functions of temperature at an applied field of 
70 kOe, where the transition occurs at 3.34 K. 



temperature phase transition may be achieved either by 
the application of a modest magnetic field of 54 kOc 3 
or a hydrostatic pressure of 1.1 kbar. 12,13 Here we have 
considered the case where the transition is approached 
by applying a magnetic field at ambient pressure. The 
field removes the degeneracy of the S = 1 triplet states of 
the dimers, and the collective singlet-triplet excitations 
separate into a longitudinal, z-polarized wave and two 
transverse modes. The collective transverse modes are 
linear combinations of propagating modes due to transi- 
tions between the ground state and the lowest and the 
highest excited states of the single dimers. In the para- 
magnetic phase, the two transverse modes are subject to 
the same renormalization effects, because the rigid energy 



shift of the excitations, 



E q = 2h, does not influ- 



ence the quantum fluctuations and n 13 = at T — 0. At 
non-zero temperatures, the thermal fluctuations imply 
that n 13 is non-zero, but the renormalized exchange in- 
teraction J xy (q) is still the same for the upper and lower 
transverse modes in the paramagnetic phase (within the 
present approximation scheme). 

The RPA renormalization parameters are determined 
to be rather substantial in the zero-temperature limit 
at zero field, because the system is close to the criti- 
cal point. The occupation number of the dimer ground 
state is calculated to be n[j = 0.935 and a = 0.393 
meV. The constant reduction of the exchange interac- 
tion is (1 + r] xy )a Q = 0.472 meV, which is about 17% 
of the maximum value of the effective exchange inter- 
action J cS (Q) = 2.8 meV. The effective singlet-triplet 
splitting is about A + a = 5.26 meV in the temperature 
range of the maximum in the susceptibility. As indicated 
by the comparison in Fig. [TJ this is in good agreement 
with that derived from the experimental data. 3,4 Hence, 
the theory is able to account for the difference found ex- 



perimentally between the (effective) energy gap of 5.68 
meV, derived from the excitation spectrum in the T = 
limit and the smaller value of the gap determined from 
the susceptibility measurements. At zero temperature 
the renormalization parameters are independent of the 
field as long as it stays smaller than the critical one. At 
non-zero temperature the lower branch E~ is much more 
easily populated than the corresponding MF level. For 
comparison, the bulk magnetization at 53 kOe predicted 
by the corresponding MF model is about a factor of 50 
times smaller at 6 K. At a constant non-zero temperature 
E~ decreases, and the thermal population of the E~- 
branch increases, when the field approaches the critical 
one. This also implies that the renormalization param- 
eter a xy increases with increasing field. When the field 
becomes larger than the critical one, a xy is reduced as 
the field is further increased. This reduction is so large 
that it is able to stabilize the ordered moment also at 
fields smaller than the critical one, i.e. the transition is 
so strongly modified that it becomes a first order one. 
This behavior of the renormalization parameters at the 
phase transition is illustrated in Fig. [5J The unphysical 
enhancement of the renormalization effects near a phase 
transition is a quite general feature of the self-consistent 
version of RPA. In the close neighborhood of the criti- 
cal point the renormalization effects show a pronounced 
sensitivity to small modifications of the model, and the 
rather good agreement between theory and experiments 
obtained here for the critical field and for the parallel 
magnetization, see Figs. [5] and [31 is somewhat fortuitous. 
Even the minor change introduced by using the exchange 
parameters derived by Oosawa et al& rather than those 



defined by Eq. 



leads to a relative increase of i 



by 



about 10% and to corresponding changes of the magne- 
tization curves and the critical field (see Fig. [5]). 

Focusing our attention on the zero-temperature limit, 
then n 13 is zero when the field is smaller than the critical 
one, but becomes non-zero in the ordered phase. As long 
as the ordered moment is small, the contributions due 
to n 13 may be neglected and the equations determining 



(S x ) and (S z 



Eqs. gO]) and gU, predict 
h 



(S 



A(noi + no3) 



(S x ) 



(43) 



when terms of the order of 9 4 and a 2 9 2 are omitted. The 
equation is only weakly influenced by the renormalization 
effects since n 01 +n Q3 ss 2 at T = 0. Nevertheless, the ex- 
perimental low-temperature results are far from obeying 
this relationship, which circumstance makes it difficult 
to reproduce the field dependencies of the two magneti- 
zation components simultaneously, as illustrated by Fig. 

m 

In their modelling of the excitations by bosons, Nikuni 



et oi. s find that (S z 



j{S x ) at T = [in this expres- 



sion we have neglected the small difference n, between 
n = (S z ) and n c , derived by Nikuni et a/., which approx- 
imation corresponds to a replacement of nrji + no3 in Eq. 
fl4"3"]) by 2]. Hence, the theory of Nikuni et al. does not 
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include the factor h/A appearing in our relation between 
(S z ) and (S x ) 2 . If this factor is included, the results for 
m xy of Nikuni et al., i.e. m ± in their Fig. 4, should be 
multiplied by y/A/h ~ 2.6 leading to a slope of m xy with 
respect to field, which is nearly twice the one derived by 
the present theory, i.e. a factor of 3-4 larger than the ex- 
perimental one (see Fig. [SJ. In the paramagnetic phase, 
the lowest excited state of a single dimer is the |1 + 1) 
state with S z — 1, and Nikuni et al. are assuming that 
this state is the one determining the wave functions of the 
lowest lying mode of collective excitations, and hence the 
one which defines the condensate in the ordered phase. 
This corresponds to assuming a — 7r/4 in our Eq. (|3T[) . 
However, as also stressed by Matsumoto et aL,— it is 
crucial to include the presence of the |1 — 1) level in or- 
der to get a consistent description of the excitations and 
of the condensate. In the paramagnetic phase, the ma- 
trix element of S x between the ground state 1 00) and the 
lowest excited state |1 + 1) is numerically the same as 
its matrix element between |00) and |1 — 1). The same 
applies to S y , and this means that the collective trans- 
verse excitations transmitted via these two operators are 
mixed S z = ±1 excitations. Being proportional to J(q), 
the degree of mixing depends on the wave vector and is 
at its maximum at the ordering wave vector. In corre- 
spondence to this, the MF ground state in the ordered 
phase, Eq. (|3"T|) . involves |1 — 1) as well as |1 + 1). The 
two states are of equal importance in the limit of zero 
field, and the relative weight of the two states is shifted 
from 1 in the presence of a field as described by the angle 
a ~ h/A. The two components depend differently on a, 
e.g. (S z ) = whereas (S x ) has its maximum at a = 0, 
and the factor h/A in Eq. (|43[) is a simple consequence 
of this difference. 

The present self-consistent RPA theory accounts rea- 
sonably well for the paramagnetic properties of the dimer 
system. Within the MF model, the bulk susceptibil- 
ity vanishes exponentially in the zero-temperature limit, 
whereas the present RPA model predicts a power law 
mjH cx with tf> = 1.8 at H = 53 kOe. This is consis- 
tent with experiments and the RPA theory also predicts 
the right critical field for the phase transition. In contrast 
to the boson model of Eq. (JTJ) , the present theory does not 
rely on any free parameters. Note, however, that if we 
use the effective exchange parameters, which Oosawa et 
al. determined from their measurements of the dimer ex- 
citation spectrum^ the critical field is derived to increase 
slightly faster with temperature than observed, as shown 
in Fig. [21 This minor discrepancy was neutralized by a 
small adjustment of the density of states of J cff (<7). Ac- 
tually, it would have been a surprise if the self-consistent 
RPA theory had been able to predict the right critical 
field without any adjustments. In all circumstances, it is 
clear that the theory needs to be corrected due to crit- 
ical fluctuations, since the self-consistent RPA predicts 
the phase transition to the antiferromagnetic phase to 
be of first order in contradiction with experiment. 

The dimer system has a number of unusual magnetic 



properties. The most outstanding one is that the system 
is driven into the phase of antiferromagnetic order by the 
application of a uniform field. Another unusual property 
shown by the ordered phase is that the bulk magneti- 
zation increases, when the temperature is lowered at a 
constant field, as this happens in spite of the fact that 
the "degrees of freedom" are being reduced because of 
the accompanying enhancement of the antiferromagnetic 
order parameter. This behavior is not in accordance with 
the MF model, whereas the self-consistent RPA theory 
accounts, at least, qualitatively for this observation. The 
self-consistent theory of the ordered phase is complicated, 
and we have been forced to neglect a number of effects. 
One of the complications, which has not been mentioned 
above, is that the matrix elements of S i between the 
ground state and the excited states are no longer zero in 
the ordered state implying additional modifications of all 
normal modes of the system (the two classes of modes 
may mix because the Cu sites lack inversion symmetry) . 
We may also add that the non-zero values of the diago- 
nal elements of S x in the ordered phase effectively give 
rise to additional contributions to A 1 and A 3 . Fortu- 
nately, these extra complications should be unimportant 
within the regime of low temperatures and small order 
parameter, where the theory is applied. 

The approximations made for b xy (u>) in Eq. ([28)) are 
only acceptable, when the renormalization effects due to 
this term are small. This term is handled in a more rigor- 
ous way by a diagrammatic high-density expansion^ To 
first order in X/z (where z is the number of interacting 
neighbors), the diagrammatic theory may be formulated 
in a self-consistent, "effective medium" fashion, which is 
the equivalent of the present, self-consistently improved 
version of RPA. 17 ' 29 This 1/z theory has been applied 
to the Ising systems H0F3 and LiHoF ^ 30 ! 31 and in these 
cases the phase transitions are predicted to remain of 
second order. Therefore, we expect that the use of the 
1/z expansion theory, to first order in 1/z, is going to re- 
produce the self-consistent RPA results derived here for 
the paramagnetic phase, and should lead to an improved 
description of the ordered phase. 

The present RPA theory establishes a classification of 
the renormalization effects, which affect the properties of 
this quantum-critical dimer system. The most important 
one is the constant reduction of the exchange interaction 
by (1 + r] xy )a xy . This term is equivalent to that produced 
by an on-site repulsive interaction J xy (ii) and is here 
found to be determined directly from the higher-order 
modifications of the RPA Green functions. Our analy- 
sis of the spin model also predicts the presence of other 
renormalization effects, such as the w-dependent correc- 
tion b xy (ui) to J xy (q) and the increase of the effective 
splitting between the single-dimer energy levels by a xy or 
a z . The q-independent reduction of the exchange inter- 
action has its parallel in the phenomenological repulsive 
interaction, v in Eq. ([T]) , in the boson model of Nikuni et 
al.^ whereas the two other renormalization effects have 
no counterpart in their theory. A more problematic sim- 
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plification made by Nikuni et al. is their assumption that 
the low temperature properties of the system are dom- 
inated by one type of bosons, whereas, in reality, the 
system contains three different kinds, where those corre- 
sponding to S z = 1 and S z = — 1 are mixed. The degree 
of mixing depends on wave vector and on field, and is 
important for the characterization of the bosons in the 
condensate. 

One basic difficulty in the many-body theory of local- 
ized spin systems is that the operators describing the dy- 
namics of the single spins are not bosonic but more com- 
plicated operators, as indicated by Eq. (TTTj) . This com- 
plication is responsible for the need to renormalize the 
simple RPA theory. Although the present self-consistent 
theory includes the leading-order renormalization effects, 



the comparison between theory and experiments within 
the ordered phase of TICUCI3 is not satisfactory. The 
diagrammatic 1/z theory , 17 ' 28 ! 29 represents a more sys- 
tematic approach and should be able to give a more ac- 
ceptable description of the ordered phase. We expect, 
however, that the pronounced experimental violation of 
the MF/RPA relation Eq. 1(35]). between the bulk magne- 
tization and the ordered antiferromagnetic moment, will 
remain a challenge to future theory. 



Acknowledgments 

We thank Kim Lcfmann for stimulating discussions. 



1 C. J. Pethick and H. Smith, Bose-Einstein Condensation 
in Dilute Gases, Second Edition (Cambridge University 
Press, Cambridge, 2008). 

2 T. Giamarchi, Ch. Riiegg, and O. Tchernyshyov, Nature 
Phys. 4, 198 (2008). 

3 A. Oosawa, M. Ishii, and H. Tanaka, J. Phys.: Condens. 
Matter 11, 265 (1999). 

4 R. Dell'Amore, A. Schilling, and K. Kramer, Phys. Rev. B 
78, 224403 (2008). 

5 N. Cavadini, G. Heigold, W. Henggeler, A. Furrer, H.-U. 
Giidel, K. Kramer, and H. Mutka, Phys. Rev. B 63, 172414 
(2001). 

6 A. Oosawa, T. Kato, H. Tanaka, K. Kakurai, M. Miiller, 
and H.-J. Mikeska, Phys. Rev. B 65, 094426 (2002). 

7 A. Oosawa, H. Tanaka, T. Takamasu, H. Abe, N. Tsujii, 
and G. Kido, Physica B 294-295, 34 (2001). 

8 T. Nikuni, M. Oshikawa, A. Oosawa, and H. Tanaka, Phys. 
Rev. Lett. 84, 5868 (2000). 

9 G. Misguich and M. Oshikawa, J. Phys. Soc. Jpn. 73, 3429 
(2004). 

J. Sirker, A. Weisse, and O. P. Sushkov, J. Phys. Soc. Jpn. 
74 Suppl., 129 (2005). 

1 M. Matsumoto, B. Normand, T. M. Rice, and M. Sigrist, 
Phys. Rev. B 69, 054423 (2004); Phys. Rev. Lett. 89, 
077203 (2002). 

2 Ch. Riiegg, A. Furrer, D. Sheptyakov, Th. Strassle, K. W. 
Kramer, H.-U. Giidel, and L. Melesi, Phys. Rev. Lett. 93, 
257201 (2004). 

3 Ch. Riiegg, B. Normand, M. Matsumoto, A. Furrer, D. F. 
McMorrow, K. W. Kramer, H.-U. Giidel, S. N. Gvasaliya, 
H. Mutka, and M. Boehm, Phys. Rev. Lett. 100, 205701 
(2008). 

4 H. Tanaka, A. Oosawa, T. Kato, H. Uekusa, Y. Ohashi, K. 
Kakurai, and A. Hoser, J. Phys. Soc. Jpn. 70, 939 (2001). 

5 S. B. Haley and P. Erdos, Phys. Rev. B 5, 1106 (1972). 

6 P. Bak, Thesis, Ris0 Report No. 312 (Ris0, Denmark). 



J. Jensen and A. R. Mackintosh, Rare Earth Magnetism: 
Structures and Excitations (Clarendon Press, Oxford, 
1991); |http://www.nbi.ku.dk/page40667.htm| 
D. N. Zubarev, Usp. Fiz. Nauk 71, 71 (1960) [Sov. Phys.- 
Usp. 3, 320 (I960)]. 

J. Jensen, J. Phys. C: Solid State Phys. 15, 2403 (1982). 
B. Leuenberger and H. U. Giidel, J. Phys. C: Solid State 
Phys. 18, 1909 (1985). 

N. Cavadini, Ch. Riiegg, W. Henggeler, A. Furrer, H.-U. 
Giidel, K. Kramer, and H. Mutka, Eur. Phys. J. B 18, 565 
(2000). 

V. N. Glazkov, A. I. Smirnov, H. Tanaka, and A. Oosawa, 
Phys. Rev. B 69, 184410 (2004); A. K. Kolezhuk, V. N. 
Glazkov, H. Tanaka, and A. Oosawa, Phys. Rev. B 70, 
020403(R) (2004). 

A. Oosawa, H. Aruga Katori, and H. Tanaka, Phys. Rev. 
B 63, 134416 (2001). 

Y. Shindo and H. Tanaka, J. Phys. Soc. Jpn. 73, 2642 
(2004). 

Ch. Riiegg, N. Cavadini, A. Furrer, K. Kramer, H. U. 
Giidel, P. Vorderwisch, and H. Mutka, Appl. Phys. A 74 
[Suppl.], S840 (2002). 

Ch. Riiegg, N. Cavadini, A. Furrer, H.-U. Giidel, K. 
Kramer, H. Mutka, A. Wilders, K. Habicht, and P. Vorder- 
wisch, Nature 423, 62 (2003). 

K. Tatani, K. Kindo, A. Oosawa, and H. Tanaka, (unpub- 
lished) . 

R. B. Stinchcombe, J. Phys. C: Solid State Phys. 6, 2459 
(1973); 6, 2484 (1973). 

J. Jensen, J. Phys. C: Solid State Phys. 17, 5367 (1984). 
J. Jensen, Phys. Rev. B 49, 11833 (1994). 
H. M. R0nnow, J. Jensen, R. Parthasarathy, G. Aeppli, T. 
F. Rosenbaum, D. F. McMorrow, and C. Kraemer, Phys. 
Rev. B 75, 054426 (2007). 



